EPJ manuscript No. 

(will be inserted by the editor) 



Perturbation theory for large Stokes number particles 
in random velocity fields 

Piero OUa^ and M. Raffaella Vuolo^ 

^ ISAC-CNR and INFN, Sez. Cagliari, 1-09042 Monserrato, Italy 

^ Dipartimento di Fisica and INFN, Universita di Cagliari, 1-09042 Monserrato, Italy. 
Received: date / Revised version: date 

Abstract. We derive a perturbative approach to study, in the large inertia limit, the 
dynamics of solid particles in a smooth, incompressible and finite-time correlated 
random velocity field. We carry on an expansion in powers of the inverse square root 
of the Stokes number, defined as the ratio of the relaxation time for the particle 
velocities and the correlation time of the velocity field. We describe in this limit 
the residual concentration fluctuations of the particle suspension, and determine 
the contribution to the collision velocity statistics produced by clustering. For both 
concentration fluctuations and collision velocities, we analyze the differences with the 
compressible one-dimensional case. 

PACS. 47.55.Kf Particle-laden flows - 46.65.-|-g Random phenomena and media 



1 Introduction 

An important component in the transport of 
aerosols by turbulent flows is the tendency to 
form clusters, an ubiquitous phenomenon that 
has been observed e.g. inside clouds [1], and 
whose simplest instance is particle aggregation 
in a ID (one dimensional) random force field [2]. 
In ID, clustering is the combined result of par- 
ticle slipping with respect to the random field, 
and the fact that the random forces pushing par- 
ticles apart become smaller as the separation 
decreases [2,3]. In more than ID, the picture 
is more complicated, and the process is accom- 
panied by preferential concentration of the par- 
ticles (supposed denser than the fluid) in the 
strain regions of the flow [4]. Thus, contrary to 
the intuition that would suggest a mixing be- 
havior, a spatially homogeneous random veloc- 
ity field will lead to the formation of clumps 
out of an initially uniform distribution [5] . This 
unmixing will take the form of concentration 
fluctuations, superimposed to a mean concen- 
tration field that remains uniform, and exeeding 
those due to discreteness in the distribution, de- 



scribed by Poisson statistics [6] . Additional seg- 
regation will be produced by inhomogeneity of 
the turbulent flow, such as in wall turbulence, 
resulting in a non-uniform mean concentration 
profile [7] . Gravity will affect the concentration 
fluctuations of heavier particles [8,9], decreasing 
the effective correlation time of the turbulent 
fluctuations sampled by the particles [10]. 

A practical motivation for the interest in 
clustering is clearly the possibility of enhanced 
binary collision, compared to a spatially homo- 
geneous condition, and this could find applica- 
tion e.g. to rain formation [11,12]. It must be 
mentioned that there have been suggestions that 
clustering is of secondary importance in the col- 
lision dynamics [14, 15]. It is not clear, in gen- 
eral, how clustering affects the relative velocity 
dynamics, and in this way also collisions [13]. 

Several models for the clustering dynamics 
have been proposed [5,16,17], however, although 
high Reynolds number turbulence is a multi- 
scale flow, the bulk of the studies has been on 
the dynamics of inertial particles in smooth ran- 
dom fields (see e.g. [18, 19]). There are good 
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reasons for this. The dynamics of a sufficiently 
small (and sufficiently dense) spherical particle 
can be described by the relaxation time of its 
velocity relative to the fluid: the Stokes time 
Ts = 2/9 rpA/fo, where ro is the particle ra- 
dius, A is the ratio of the particle to fluid den- 
sity and is the kinematic viscosity of the fluid 
[20]. Now, most atmospheric aerosols arc char- 
acterized by Stokes times shorter than the Kol- 
mogorov times of the flows by which they are ad- 
vccted [11]. Furthermore, experimental data [21] 
and numerical simulations [4] both indicate that 
clustering is stronger for particles with rg of the 
order of the Kolmogorov time, corresponding to 
a range of scales in which the velocity field could 
be approximated as smooth. 

In the case of a smooth incompressible ran- 
dom field, clustering is peaked at S = ts/te ^ 
1, where te is the correlation time of the field 
and S is called the Stokes number. For an in- 
compressible velocity field, the limit S ^ cor- 
responds to passive scalar transport, and spa- 
tially homogeneous initial conditions will not 
lead to clustering. In the opposite limit S oo, 
one expects that the particles be scattered by 
the velocity fluctuations they cross in their mo- 
tion [22] as if undergoing Brownian diffusion, 
resulting again in no clustering. One reason for 
being interested in this limit is that also parti- 
cles in turbulent flows for which rg lies in the 
inertial range (or above), see smaller vortices 
as if 5* ^ 1. It must be said that the large S 
limit is more of interest for industrial applica- 
tion, since in cloud turbulence, the relative mo- 
tion of heavier droplets is dominated by the dif- 
ferent gravitational settling velocity of particles 
of different size [11] and the possibility of chaotic 
trapping [23,24]. 

We shall focus in this paper on the case of a 
strongly stirred fluid, such that, also for large-S* 
particles, the largest component of the particle- 
fluid relative motion is due to inertia and not 
gravity. We will try to understand, in particular, 
how the Brownian limit of [22] is achieved. 

There is some evidence [25] that clustering 
destruction occurs at finite S, as the result of 
a crossover of the correlation dimension D2 of 
the particle distribution, above the dimension 
of the space D. For S below this threshold, the 
probability density fmiction (PDF) p{r) of find- 
ing a pair of particles at separation r diverges 
for r ^ like r^^-D (gee also [18,26]), while 
it should remain finite above thcshold. Actu- 
ally, residual (non-singular) concentration fluc- 



tuations remain present above threshold, decay- 
ing in ID like p{r = 0) oc S~^^^, and their origin 
lies in the fact that only particles moving at in- 
creasingly small relative velocities, as 5* ^ 00, 
are able to stay close long enough to remain cor- 
related [13]. 

The picture just described suggests the pos- 
sibility of a pcrturbativc expansion around S 
00. This is formally identical to an expansion 
in powers of the Peclet number Pe, for parti- 
cle advcction in the presence of strong molecu- 
lar diffusion. In both the large 5* and the small 
Pe regimes, the small quantity is the contribu- 
tion to the relative particle motion produced at 
small separation, by the correlations in the ran- 
dom field. The perturbative approach will allow 
to carry on an analysis of the large-S" inertial 
particle dynamics for D > 1, where the heuristic 
approaches like the one in [13] are made difficult 
by the complexity of geometry and incompress- 
ibility. 

This paper is organized as follows. In Sec. 
2, the model equations for the problem are de- 
rived. The perturbative approach is derived in 
Sec. 3. In Sec. 4, the approach is tested in ID, 
comparing with numerical simulations and with 
the heuristic results in [13]. In Sec. 5, concen- 
tration fluctuations are analyzed in the case of 
a 3D incompressible random field. Section 6 is 
devoted to determination of the effect of con- 
centration fluctuations on the relative particle 
velocity statistics. Section 7 contains the conclu- 
sions. Discussion of additional technical aspects 
is left to the appendices. 



2 The stochastic model 

Following the derivation in [13], we introduce 
a zero mean, smooth Gaussian random velocity 
fleld u(x, t), with correlation 

(«„(r, t)u0{0, 0)) = alF{t)gc,0{r), (1) 

where gafjiO) = Sa/j, F{0) = 1, and we can in- 
troduce correlation times and lengths te and 
obeying /°° dr F(r) = te, and drgaffirr) ~ 
Tv (f = r/r fixed). We assume isotropy and ho- 
mogeneity of the field in space and time. 

From the field parameters Tv, te and (T„, we 
can introduce the Kubo number K = auTE/r^, 
which tells us whether the field is intrinsically 
short- or long-time correlated. The K ^ limit 
would correspond to a random field with zero 
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correlation time, like in the Kraichnan model 
[27]. The K ^ oo limit would correspond to a 
frozen field regime. In most of our analysis we 
shall assume K — 0{1), as in actual turbulent 
flows. 

A suspension of inertial particles is advected 
by the random field and their velocity v is taken 
to obey the Stokes equation: 

V = Tg M-v(i) + u(x, i)] + 77, x = v, (2) 

where we have allowed for a Brownian motion 
component rj: {ria{t)rii3{0)) = K,Tg'^5ap5{t)] k is 
the molecular diffusivity of the particles. 

We choose units such that (t„ — ts = 1; 
therefore: 

TE = S-\ rv = {KS)-\ 

In the regime S ^ 1, the particle displacement 
in a correlation time te becomes negligible with 
respect to and Eq. ([2]) could be approximated 
by a Langcvin equation (the precise condition 
for a Langevin dynamics is ^ S [13]). It is 
then possible to substitute into Eq. 

u(x(t),<)^(2rB)i/2|(t), 

with ^{t) white noise: {£,a{t)£,i3{0)) = SapSlt), 
i^aVls) = 0. We thus obtain for the particle ve- 
locity variance ~ S~^{1 + Pe~^), where the 
Peclet number Pe = 2a'^TE / k parameterizes the 
relative strength of the random field and molec- 
ular contributions to particle diffusion. 

Turning to the relative motion of particle 
pairs, let us introduce difference variables 

IV = V2-Vi, r = X2-Xi, 

where 1 and 2 label members of a particle pair. 
(The particles are assumed immaterial, so that 
they can cross without interaction). As in the 
one-particle case, we can approximate the equa- 
tion for the relative motion of particles with the 
Langevin equation: 

i'a = -Va + baf3{r)^f3, Ta^^a- (3) 

where 

bai{'c)h^l3{r) = q^[5ai3 - agapi^)] (4) 

(in this paper, summation over repeated indices 
is assumed), with 

a= [1 + Pe)-'^Pe and q = 2{Sa)-'^/^ . (5) 



We see that a goes to zero with Pe 0, and, 
as expected, the relative particle dynamics de- 
scribed in Eqs. (I3]|4]) becomes uncorrelated. The 
zero molecular diffusion limit corresponds in- 
stead to Pe — !■ oo and a = 1. 

It should be mentioned that, in the Langevin 
equation limit, the dynamics becomes equiva- 
lent to that of a Kraichnan model and can be 
described in terms of the single parameter [3] 

e = K^S. (6) 

To understand the connection with the Kraich- 
nan model, notice that e^^Ts ~ rl/{a^TE) is the 
diffusion time of a single passive tracer across a 
distance and plays the role of effective cor- 
relation time of the field, so that e becomes an 
effective Stokes number. More interestingly, for 
Pe,e S> 1, e~^^^Ts is the permanence time at 
separation < — {KS)~^ of a particle pair 
approaching at relative velocity (T„ ^ S~^/'^ (it 
is easy to see that such particles cross at large 
enough speed to behave ballistically at scale r^). 

The limit 5 ^ iiT ~ 1, beyond allowing a 
Langevin equation based model, leads to de- 
coupling of the difference variables {v, r) from 
the center of mass variables ^(xi -I- X2) and 
i(vi -l- V2). At equilibrium, this means that the 
two-particle PDF will be in the form, using spa- 
tial homogeneity: 

p(vi,2,xi,2) = i7"V(-^^^-y^)p(i^,r), 

where >!7~^, with Q the volume of the domain 
for X, is just p{[xx+X2)/'2,). (Where not ambigu- 
ous, we do not use subscripts to identify PDF's 
of different quantities). Multiplying by iV^, with 
N the total number of particles in the domain 
and integrating over d'^uid'^W2, we obtain the ex- 
pression for the concentration correlation: 

(n(r)n(O)) ^n'^Q j d^vp{v,Y), (7) 

where fi = N/ fl is the mean concentration. The 
quantity /(r) = i7p(r) — 1, with/9(r) = / i\^vp{u, 
r), gives the strength of the concentration fluc- 
tuations (we have exploited isotropy). The con- 
centration variance can be expressed in terms 
of the function / by means of the relation ((n — 



4 



OUa and Vuolo: Perturbation theory for large S particles in random fields 



3 Perturbation theory 

For small Pe, the two-particle PDF p{i',Y\t) = 
{6{u{t) — v)5(v{£) — r)) can be determined solv- 
ing Eqs. ([SlU]) for v(t) and x(t) pertm'batively in 
a. We shall argue and verify numerically in the 
next section, that the same perturbative strat- 
egy remains valid for 5* ^ if ^ 1 also in the 
absence of molecular diffusion, i.e. for a = 1. 

The perturbative expansion of v(t) and v(t) 
is obtained substituting the Taylor expansion 

-^5a7(r(i))57/3(rW) + ■■•], 

which is obtained from Eq. ^ , into the solution 
ofEq. ©: 



-T-t 



+ /!^dr(l-e-*)6„^(r(r)K^(r), 

and Taylor expanding again in a. The ground 
state of the expansion is therefore: 



L(0) 



-T-t 



) (8) 



+g/!^dT(l-e--*Ka(T), 
while the first two corrections are: 



^SW = -(W2)5a/3(r<°'W), 
and 

b'^^^it) = -q[{a/2y;Ht)d,g^f,{r^->{t)) 

-)-(aV8)5o^(r<")(i))57Mr<'"W)], (10) 
^L^'W =/!Tdr(l-e^"*)6i^^(r)fM^), 

where = d/dr^. We focus here on the PDF 
for the particle separation r, leaving the analysis 
of the joint PDF p(iy, r) to Sec. [i The PDF for 
r, given some distribution of initial conditions 
(r(— T), i/(— T)), will be something in the form: 

p(r;i) = / ndV'-V(K'};i)<5(E'-''' -'•)■ 

•' k k 

Depending on the choice of initial conditions for 
the trajectories ending at r, p(r; t) will be in 
general a non-equilibrium PDF. Carrying out 
the integral over r'"' : 

p{v- 1)^ JY[ d^f (''V({f r<"' = r - f ; t), 



where f*''' — r<'°' for k > and f = X]fc>o 
Taylor expanding in a leads to a perturbative 
series p = p'"' + p'^' + . . ., with p^"'{r,t) = 
p{r'-"\t) = r). In explicit form: 

p(r,i)-[l-a„(r<,i)|r)-a„(dr) 

+ (l/2)9„90(Cr<^'|r)-f ...]p<°)(r,t) 

where (r'^'lr) is a shorthand for 



'(-T) = 0) (11) 



and similarly for the other conditional averages. 
Sending T ^ oo and using limy^cx) p(r*°') = 
p<")(r) = f2~^, we obtain the perturbative ex- 
pansion for the equilibrium PDF: 



f2p{r) = 1 -a„(r(,^'|r) -9„(r<,^>|r) 



+ (l/2)a,a^(r<^)r<^'l 



(12) 



Clearly, the dependence on the condition r<^>( 
— T) = in Eq. (|11|) disappears when T — )■ oo, 
but the limit must be taken after the average. 

In order to take into account the simultane- 
ous conditioning in the past and in the present 
in Eq. (fTT|) . the following expression for (r'^'jr) 
can be utilized: 



'(-oo) = 0), (13) 



where the factor f2 is simply [p(r'°' (i)|r*^' (— oo) 
= 0)]-^ Substituting Eqs. dH) and ^ into ^ 
and setting without lack of generality i = 0, we 
obtain: 



dr(l - e-)(6L^> (r)eMr)'5(r''"(0) - r)). 



Using the functional integration by part for- 
mula [28], allows to treat the correlation be- 
tween £,i3{t) and the other factors in the inte- 
gral: 

^^dr(l-e-)(&l^^(r)^^9,o,(„)<5(r-'(0)-r)), 
with 6/S^p{t) indicating functional derivative; 



from Eq. 



dr^S\0)/6^piT) = <z(l - e^S,}, 



[notice that (5r<"> (t)/(5^,3 (r) = 0]. Substituting 



(0) 



5 (0) 

nally: 

/^(i 



-d^ and using Eq. we obtain fi- 



\r) = qq'ad0j°^dTil-en' 
x(5„;5(r('"(T))J(r(°'(0)-r)). 



(14) 



Similar expressions can be derived for the higher 
orders in the expansion for p, and they will be 
required to take into account incompressibility 
in the D > 1 case. 
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4 Small Peclet vs. large Stokes 

The two-particle distribution in Eq. is ex- 
pressed as a formal (regular) perturbation ex- 
pansion in a. The real expansion parameter, how- 
ever, is the ratio r/rv, where f ~ r*^' is the con- 
tribution of the velocity correlation g^is to the 
relative particle displacement. A relevant quan- 
tity is then the permanence time at r'°' < r^, 
i.e. the time during which gap{r^"^) is signifi- 
cantly different from zero and can contribute to 
the integral in Eq. (jl4p . If r < r^, the perma- 
nence time will be r-v/i^ ^ r^/q [see Eq. ([8])]. 
Taking d^^dp ~ r^^ and S{r<°^ - r) - f2-\ 
gives therefore: 

r/r - a„(ri^'|r) ^ (aK/q) min(l, (r^/q)^). 

Now, our perturbation expansion is carried on 
in powers of a for fixed q; in other words, a 
(and therefore Pe) are considered small inde- 
pendently of q. Had K and 5, in place of q, 
been chosen as the fixed quantities in our prob- 
lem, we would have obtained instead, from Eq. 
®: 

f/rv ~ a^/^e^/^ niin(l, a/e). (15) 

This suggests that the perturbation expansion 
of Eq. (fT^ could be translated, for large e and 
Pe generic, into a new expansion in powers of 
e-V2. 

The validity of the new expansion rests on 
regularity of the solution p(r) in r = (the es- 
timate da, dp ^ adopted imphcitly this as- 
sumption). While this is guaranteed for small a 
by molecular diffusion, for a = 1, baf3{0) = 
and Eq. ^ becomes singular. Now, the corre- 
spondence between the small a and large e cases 
guarantees finiteness of the terms in the series 
(it is possible to see that this is not verified in 
the case a ground state of ballistic particles at 
scale Tv is adopted [13]). Hence, expanding in 
powers of e~^^'^ is a meaningful procedure. It is 
still possible that the perturbative approach is 
simply unable to catch singular behaviors [no- 
tice that Eq. does not reveal break-down of 
the perturbation theory for a = 1 and e ^ 1, 
in spite of the singularity of p(r)]. The numer- 
ical evidence in [25], that singular behavior are 
absent for large e, however, suggests that this is 
not the case. 

Let us test in ID the perturbative approach 
just introduced. For the sake of definiteness, we 
consider a Gaussian profile for {u{r,t)u{0,t)): 

g{r) = exp(-rV2r2); (16) 



hence: b'"' = q and 6<^' = {qa/2) ^(rC"). 

Substituting Eq. p6| into (fT4|) and then into 
()12p . we can write p'^'(r) in the compact form: 



{t)Ar{r,t), (17) 



with Ar{r,t) = fdr^°>{t)p{r^°\t),r^'»{0)^r)g{ 
r'-°''{t)) and Hr{t) = (1 — e*)^(}^; the averages 
have eliminated all dependence on the initial 
conditions att — — T. We can write p(r'''> (i), r<"> ( 
0)) = l?-V(r<'"(*)k*°'(0)) = n-'^p{s{t)), with 
s{t) = r""(i) - r<">(0), and p{s{t)), for T ^ oo, 
is a zero mean Gaussian with variance, from Eq. 



's{t) 



q^[~l-t + e% 



(18) 



Substituting into Eq. ^T7} . we obtain the result 



fir) 



4P72 Jo 
1 



G.(e,t) 



exp 



{ 2g7(m)}' 



(19) 



Grie ,t)=t-l 



where r = r/q. Notice that q is the typical sep- 
aration of a pair of particles at a time ^ rs = 1 
after crossing (f played the role of outer scale 
in the matching asymptotics analysis in [13]). 
In other words, the correlation length of the 
concentration fluctuations is not expected to be 
of the order of r^, rather, of the larger scale 
qr^{e/ay/^r^. 



4.1 The ID case: small e regime 

Taking the limit e ^ under the condition 
K^/S <C 1, guarantees that the Langevin ap- 
proach of Sec. [2] continues to be valid. (Notice 
that the above conditions imply K <^ 1). This 
regime has been extensively studied in [3, 19] as 
a stochastic model for the behavior of low iner- 
tia S ^ 1 particles in realistic turbulent flows, 
for which K ^ 1. 

The e limit is diffusive at scale for the 
variable r{t) and we can derive an exact expres- 
sion for p{r), valid for all a < 1. This expression 
can then be used to test the perturbative results 
of the previous section. A diffusive limit means 
that the variable r{t) changes in a correlation 
time Ts much less than the spatial scale [on 
the contrary, for most pairs, the large e limit is 
ballistic at scale and becomes diffusive only 
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at the scale (e/a)^/^rv]. The drift and diffusiv- 
ity in the effective equation for r are the leading 
contributions to 

([r(i)-r(0)]|r(0))A and ([r(t) - r(0)]2 |r(0))A, 

for t > 1 and \r(t) - r(0)| < r^. Actually, it 
is easy to prove from moment analysis of the 
Fokker-Planck equation for Eqs. ([3]), that the 
drift (r|r) = (i^|r) is zero [13]. We only have to 
calculate the diffusivity; from Eq. we can 
write, for t 3> 1 and \r{t) — r(0)| <C rv- 

r{t) = r{0) + 6(r(0)) / drCir) + 0{t°) 
Jo 

and {[r(t) - r{0)]^\r{0)) /t ~ B{r{0)). The par- 
ticle separation behaves therefore like a non- 
uniform Brownian motion: f = b{r)^ and from 
the associated Fokker-Planck equation we ob- 
tain the exact expression: 

fir) = [l-aeM-liKSrf)]-\ (20) 

Notice the quadratic divergence of the PDF at 
r for a = 1, that is approached at i — > oo 
with an 0(e) rate [3]. 

Let us compare with Eq. (fT9|l . In the hmit 
e ^ 0, we can write Gr{e,t) = t + a/(4e) and 
Eq. (I19|) becomes, after a few manipulations: 

f{r)=aexp{-^{KSrf), 
but this is precisely the 0(a) contribution to Eq. 

mi. 



4.2 The ID case: large e regime 

The prediction in [13] that /(O) - e^^/^ for 
00 seems to be supported by Eqs. (fTSj) and 



([TO)) . However, the integral in Eq. pl?|) presents 
singularities for e ^ oo, and we must put some 
care in the analysis. 

Let us consider first the case r — 0, and write 
Eq. (fT9|) in the form 



/(0) = 



7,5/2 . 



t^dt 



(^2 + a/(2e))3/2 



where e^^/^ ^ T ^ 1; the dots indicate the 
remnant of the integral in Eq. (|19p . which is 
shown by inspection to be finite in the limit 



e — i- oo. The first integral in the formula above, 
instead, is dominated by the ballistic crossing 
time t ^ iT^I"^ and is logarithmically divergent 
for e — !■ oo. We obtain the leading order expres- 
sion for /(O): 



/(O) 



a5/2 Ine 



2V2e 



The logarithmic divergence is eliminated and a 
pure power law is recovered, provided we coarse 
grain the function / at a fixed scale i?: 

Uir) = (l/R) f w{r'/R)f{r-r')dr', 



where w{r) is a smoothing function with w{r) > 
and J w{r)dr = 1. A simple analytic expres- 
sion for fn is obtained choosing a Gaussian w(r) 
= (27r)-i/2cxp(-rV2): 



,5/2 



(1 



-t\2 



fdt 



4Vi 7o [t~i 



i?2]3/: 



r, (21) 



where R = R/q. 

These results confirm the heuristic predic- 
tion in [13]. The interesting point is the im- 
proved performance of perturbation theory at 
large e, as illustrated in Fig. [TJ The collapse of 
the rescaled profiles indicates a region in which 
the higher orders in the perturbation expansion 
are playing a negligible role. The perturbation 
expansion in a can therefore be converted into 
one in e~^^'^, as claimed. 

We focus next on the correlation profile /(r). 
We can obtain analytical expressions for large 
r = r/q. We rewrite Eq. pop in the form 

5 /2 /"OO 

f(r) = / f*h^(x)U(x + h{x))dx, 

4Ve Jq 

(22) 

where x = t/r'^, h = r 2[exp(— r^a;) — 1] and 
U{x) — x'^^"^ exp{—l/{2x)). We can Taylor ex- 
pand U{x + h{x)) = Ulx) + h{x)U'{x) + .... We 
then substitute into Eq. ([22]) and integrate by 
parts, using the fact that R and all its deriva- 
tives are zero at both zero and infinity. The in- 
tegrand in Eq. ([22]) takes the form: 



r'[h'^{h^y + ^{hr^-]U{x) 

= [1 + aic-^'" + aze-'"'" + ...]C/(a;). (23) 

We have ai = -|f + - ^ + ... ^ - exp(-l) 
and the leading order in Eq. (I23p is therefore [1— 
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10000 



Fig. 1. Rescaled fluctuation amplitude profiles for 
different values of a, coarse grained at scale R — 
0.08. Cases a—d are from from numerical integration 
of Eqs. daill: (a) a = 1.0; (6) a = 0.9; (c) a = 0.7; 
(d) a = 0.4. Case (e) is the theoretical prediction 
from Eq. ((2l|. The heavy line is e"^/^ 



exp(— 1 — f^2;)]x exp(— l/(2x)). Substituting 
into Eq. (j22|) . we get the final expression: 



fir) 



,5/2 



exp 



(-A/2|f|-l)+0(e-2|^l). (24) 



Notice that, opposite to the case of /(O), the 
integral in Eq. ([22]) is dominated by the long 
time scale t ^ f, corresponding to x f~^. The 
heuristic predictions on the correlation profile 
in [13] are therefore confirmed. As shown in Fig. 
2, however, agreement with the scaling predic- 
tion of Eq. p4|) is obtained only for rather large 
values of e and of the rescaled separation f. 



5 Concentration fluctuations in 3D 

From now on we restrict the analysis to the 
regime Pe — > oo (i.e. a = 1). This is physically 
consistent with a large S regime, correspond- 
ing to particles with high inertia in the presence 
of strong advection. We recall that for K ^ 1, 
S ^ 1 corresponds to e 3> 1. In more than 
ID, if the flow is incompressible, concentration 
fluctuations will appear only at second order in 
the perturbative expansion of Eq. (|12p . In fact, 
the generalization to more than ID of Eq. (fT7|) 
reads: 

p'^'(r) = -(2/5)/°^dt(l-eT 

x9„50/dV«)(%„/3(r<"'(t)) 
xp[r<''>(i))|r('''(0) = r]. 




0.01 



Fig. 2. Fluctuation correlation profiles for a — 0.9 
and three different values of e. Dotted line: e — 10; 
heavy line: e = 100; thin line: e = 1000. Diamonds 
correspond again to the e = 100, a = 0.9 case, but 
from direct simulation of Eq. ((2]), instead of Eqs. 
([3]|4]). Simulation parameters in the case of the full 
simulation: K = 1, N — 10* particles, total simu- 
lation time = lOOrs; domain length f2 = 105'"^''^, 
corresponding to 512 modes in the Fourier decom- 
position of the random field u{x, t). The straight line 
is the /(r) oc exp(— \/27^) prediction of Eq. (|24p . 



that will be identically zero if dagap ~ 0. The 
lowest order contribution to /(r) is therefore: 



1 



where r'^> and r'^> are given in Eqs. ([9llT0)). Con- 
trary to the heuristic prediction in [13], of an 
e~^/^ decay at large e for the concentration fluc- 
tuation variance, we thus expect an behav- 
ior. We show that this estimate is correct in Ap- 
pendix A, by explicit calculation of p*^' . 



6 Relative velocity statistics 

The approach in Sec. Ill can be extended to the 
calculation of the joint PDF p(r, i/), with Eqs. 
()12I13I) being replaced by 



p{y) = [i~d^{y^-^\y)~dM:'\y) 



idMy'^%'\y) + --¥"\y) 



(25) 



and 



'|y)p""(y) 

y<^)(t)J(y<")(t)-y)|y<^)(-T) = 0) 
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where y = (r, v) and the vector indices run now 
from 1 to 6. Clearly, integrating Eq. p5|) over u 
leads to Eq. (fT2l) . 

From Eq. ([^5]) . we could in principle obtain 
information on the velocity components that con- 
tribute the most to clustering. Looking at Eq. 
(PS)) . we see that the terms contributing to p*'"' (r) 
for fc > are those containing only spatial deriva- 
tives (the others give zero after integration over 
f). We can thus write: 

p(r, u) = [p^°\u) + Pciii^lvMv) + n-'pBii^\r), 

(26) 

where f2 ^pB{i^\r) contains the terms in the 
RHS of Eq. ((25)) involving velocities derivatives 
= d/dys+a that do not contribute to Eq. 
(fn)l : of course, / d^iy pB{i^\r) / d?v pdiu |r) = 
0. Substituting into Eq. ((^ : 



Integrating Eq. ([5]) with initial conditions [u = 
0, r = 0), we see that the particles separate at 
r<°> = Tv at a time ~ e~^^^Ts, which is the low- 



l-aa(r<,^>|y) 



<i,^a,| 



(27) 



p,M^)^[[Qp{v)]-^ 

+ •••]- l}p'"'('^)- 

We may thus interpret pci{i'\v) as the clus- 
ter contribution to the velocity PDF p(t'|r) at 
separation r. It is important to notice that this 
contribution is not necessarily localized in the 
clusters: different moments with respect to u of 
Pc/(^',r) — pci{i>\v)p{r) do not come necessarily 
from the same spatial spots in the volume J7. 

Contrary to Eq. (jl2p . the expansion param- 
eter in the velocity part of Eq. ((25)) is now ve- 
locity dependent. In analogy with Eq. ((T4)) . we 
can write 



est order estimate for Tea^it (0|0). From Eq. ((29)) . 
t^^/^Uv appears to be the relative velocity at 
t — Tea;i((0|0), and is therefore the escape veloc- 
ity out of the correlated region r < for slow 
particle pairs, that do not behave ballistically 
in that region [13]. Notice that there are no sin- 
gularities &i V = and, in contrast with a per- 
turbation theory with a ground state of ballistic 
particles, the individual terms in the expansion 
remain finite as — > 0. 

Here, we focus on the contribution to the col- 
lision velocity variance (i^^|r=0)c; = / d^vv'^pci{ 
h>\v — 0). Contrary to perturbation theory for 
p{i'), that breaks up for v ^ e~^l^a^, the one 
for the moments of p is perfectly well behaved, 
as the contribution from v < (T^I^ to (i/P|r)c; — 
J d^iy lyPpciMr) is - e-(3+p)/6^^^(o|^)^ ^j^d p,i{0 
|r) is finite. 

Substituting into Eq. (|27p. we obtain 



Pel 



y\r),i = -a„(r<,^)(z.('")2|r) -9„(C(i.(°')2|r) 
+ ia^a^(r<^'r<^'(z.<«))2|r) + ... 
-/(r)((i.<"))2)-f /(r)a„(C(i.<'")2|r) 

(30) 

and we recall that /(r) = np{r) — l [see Eq. ([7])]. 
Again, all conditional averages are intended in 
the sense of Eq. (jlip. Analogously to the anal- 
ysis in Sees. IV and V, we see that the first 
non-zero contribution in the incompressible case 
arises at second order: 



x(g%(r<''>(r))%"')(0)-y)). (28) (z^^jr)'^' = -9„(r<^'(:.<°')2|r) - /2p'^)(r)((i.<«))2) 



2^3 



(z.|r)p<°)(i.) 



+ i9^9^(r<^V(i'(:.(°))2|r), 



where rea;it(i'|r) is the permanence time at sep- 
aration < Tv of a particle pair, characterized at 
given time by values (r, v) of the relative posi- 
tion and velocity; assuming ballistic motion and 
taking r = 0: Tgxit ^ T-^jv. For typical particle 
pairs, for which j/ ~ CTi, ^ S*"^^^, we thus have 
p(i) ^ e~i/2p(o)^ Pqj. ^ _ e~i/^cr„, though, we 
find p*^' ~ p*"' and perturbation theory breaks 
down. 

To understand what happens, we derive from 
Eq. ([3]) the analog of the lowest order equation 
for the position ([5]), in the case of the velocity: 



while in the ID case: 



(31) 



(29) 



(32) 

The calculation of the cluster contribution 
to the collision velocity variance at this point is 
a matter of lengthy but straightforward algebra. 
Leaving the calculation details to Appendix B, 
the result in ID and in the incompressible 3D 
case are [see Eqs. ((B2IB6I) ]: 

((i.<">)2)-i(z.2|o)<;) = -4(^e)-i/2lnc^, (33) 

and 

((:/<°')2)-i(t.2|^ = 0)^^' =i/e-i/2, (34) 
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with c and H positive 0(1) constants. 

We may interpret these results by saying that 
concentration fluctuation produce collision hin- 
dering in the compressible ID case, and to colli- 
sion enhancement in the incompressible 3D case. 

Notice the 0{e-^/^) behavior of {v'^lr = 0)^"' 
in place of the expected 0{e~^) at second or- 
der, ft is possible to show, however, that the 
0{e~^^^) contributions from pd and pB cancel, 
so that the total correction to the collision ve- 
locity variance is 0(e^^). 



7 Conclusion 

We have derived a perturbative approach for the 
two-particle statistics of a randomly advected 
inertial particle suspension, that is valid both 
in regimes of small Peclet number Pe (strong 
molecular diffusion) , and large Stokes number S 
(corresponding to high inertia), fn both cases, 
one expands around a lowest order of indepen- 
dent Brownian particles. This is natural in the 
small Pe regime; for large S [more precisely, for 
large e, with e defined in Eq. the mechanism 
is more subtle and is due to the fact that par- 
ticle trajectories evolve on a characteristic scale 
e^/^Tv, much larger than the correlation length 
Tv of the field u. The particle pair trajectories 
are dominated therefore by large separation, un- 
correlated u contributions [22]. 

The perturbative approach we have derived 
provides an analytical description of how the 
Brownian particle limit of [22] is achieved, that 
goes beyond the qualitative considerations in 
[13]. The concentration fluctuation amplitude 
n~^((n — n)'^), with n the mean concentration, 
appears to be 0(e^^/^) for compressible, and 
0(e~^), for incompressible flows. The correla- 
tion length of the fluctuations is ~ e^/^rv, in 
contrast with the e <C 1 case, in which, a power 
law at separations r < would occur. 

The perturbative approach allows to iden- 
tify a concentration fluctuation contribution to 
the statistics for the collision velocity v, in those 
terms in the expansion for the joint PDF p{i/, r = 
0), that lead to deviations from the uniform 
fluctuation-free regime in the separation PDF 
p{r). This goes beyond the observation that the 
expected relative velocity should decrease at small 
values of the separation r. 

It has been noted in [29] that for e <C 1, 
different collision velocities originate from par- 



ticle "jumps" starting at different initial separa- 
tions. The collision velocity distribution is thus 
affected by the spatial structure of the clusters, 
which produce an effective collision hindering. 

In the present large e regime, no such direct 
physical association between clustering and col- 
lision dynamics exists. From an analysis of the 
collision velocity variance, we see that collision 
hindering occurs in ID, while enhancement oc- 
curs in 3D if the flow is incompressible. The re- 
sult in ID is not unexpected, as clusters are in 
this case the result of particles slowing down 
relative to one another as they get closer. The 
role of incompressibility and 3D in leading to 
collision enhancement is less clear. 

Additional information on the velocity statis- 
tics is contained in the velocity structure of the 
concentration fluctuations. From Eq. we 
see that, for D < 3, the integral p(r) — J d^i/pi^ 
V, r) is dominated by small relative velocities, 
while for 13 = 3 all velocities contribute equally, 
down to the velocity scale e~^^^av, at which, 
particle pairs with r < cease to behave bal- 
listically. (The Brownian regime at large e is cor- 
related in time at the scale of the Stokes time 
T5, and one expects particle pairs with a rel- 
ative velocity that is not too small, to behave 
ballistically at separation r < r^). 

At least in ID, the predictions of perturba- 
tion theory begin to be valid only for rather 
large values of e. The fact is that clustering 
at e < 1, and residual concentration fluctua- 
tions at e 3> 1, are very different in nature. 
At small e, in first approximation, the parti- 
cle separation r, for r ^ evolves as a dif- 
fusion process with diffusivity oc er^. Clustering 
arises therefore from trapping at r = of the 
particle pairs. For large e, most particle pairs 
at r < Tv behave instead ballistically. In this 
regime, concentration fiuctuations and associate 
particle velocity modifications, can be seen as 
corrections to ballistic motion at scale Tv. They 
are at most only the remnants of the trapping 
behaviors that dominate the pair dynamics at 

Appendix A. Calculation of p<^'(r). 

The generalization of Eq. to the calculation 
of (r'^'jr) and (r'^'r'^'jr) is obvious. The result- 
ing integrals in the form 

J dtidi2 . . . {f[^; hM. ■ ■ ■)i{ti)i{t2) ■■■) 
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are simplified by repeated application of the func- 
tional derivation by part formula {F[^]£,a{t)) — 
{SF[^]/S£.a{t)), with the relations S^a{t)/S£.i3{t') 
= 6^f,S{t- t') and 5r^°\t)/5ip(t') = q5^pe{t - 
t'){l — e* ~*) with 6{t) the Heaviside step func- 
tion [e{t > 0) = 1; e{t < 0) = 0]. After some 
algebra, we obtain the result 

np^-\v) - (gV4) 

x/°^dt„/!"^dt,(l-e*")V(*-*») 

xp[r^»^ita) = r,,r('"(4) = rfc|r<")(0) = r]. 

(Al) 

As in the ID case, we expect that the time in- 
tegrals be dominated by \t\, \t'\ <^ 1 {ts — 1), 
so that (1 - e*i)'*e2(*"~i) t\ and the corre- 
lation matrix entering the Gaussian joint PDF 
p[r'">(i,) = Va,r'"\h) = rtlrCHO) = r] read: 

{K\ta)-r^WP{h)-rp]) 

= {q^/2)5o,p{l-\ta-ti,\/2)tah. 

Exploiting isotropy and incompressibility, it is 
convenient to write the random velocity corre- 
lation in the form: 



gapir) = rl[dadi3 - 5ap^'^]C{r /r^ 



(A2) 



We see that C{x)5ap is the spatial correlation 
for the vector potential for the field u and the 
Fourier transform Ck = J d^x e~^^'^C{x) is thus 
positive defined. To calculate the concentration 
fluctuation variance, we set r — 0. Writing in 
terms of Fourier components and using Eq (|A2p , 
we can then rewrite the right hand side (RHS) 
of Eq. (|A1[) in the following form: 



and 



I-oottdtat^dUj 



X[(k, • k,)2 - (fc,fcfc)%"k!t(ka,kb), 



kbr^ 



Zf^4j,(ka,kf,) = exp 



^((iafca)' + (4fc6)' 
-2t,4(l-i(i„-i6))ka-kfe) 



is the generating function for r<">(ia.6) condi- 
tioned to r'"*(0) = 0. The multiple integrals in 
p'^' are simplified passing to polar coordinates: 



(e^/^ti, 6^/2^2) = (r^/^ coscp,ri/2 sin^s) and k 
k2 — kik^z. We are going to verify that the in- 
tegrals in are dominated by s, t ~ 1 so that 
the term i(ta-4))ka -kf, in Zt^tA^a-, kb) can be 



cost 



,1/2 



smf 



1 • 



disregarded. In this case, the r integral can be 
carried out explicitly, and we obtain the result 

f{o)^^j:'dsrj'dej\dzrj'A^ 

8^(1-2^)2 sin** 26* cos* ^ Cj^^Cf^^ 
[l+cos(2(e+v3))-|-(l+z) sin 26 sin 2ip]3 ' 

(A3) 

From Eqs. (|llA2p . we have that C^.<;^ ^ 1, so 
that, provided the integrals in the RHS of Eq. 
(|A3p converge, we have /(O) oc e^^. Conver- 
gence is also sufficient to verify correctness of 
the ansatz s, r ^ 1 in the integral. It is suffi- 
cient to prove convergence near the singularity 
a.t 9 = 7r/2, = 0, z = —1, where the integrand 
is oc e^z'^ +i^2-)-3 ^^^Yi ip ^ ip + ~ 7r/2, 
6 = n/2—9 and z = l+z. Integrating first in d(f, 
leads to an expression whose leading term in z 
and ^ is oc 9z~^^'^, the remaining integrals con- 
verge as required, and therefore p'^' = 0(e^^) 
as expected. 



Appendix B. 
ance 

B.l The ID case 



Collision velocity vari- 



Substituting Eqs. ^ and ^ into Eq. ([32]), we 
easily obtain: 



= ^/;dn/:dr,/:dr3(l 
)e-^+-3 (9^,0,5(1)^(1)^(2)^(3) 



(Bl) 



X(5(r('"(0) - r) 



2^3 



where g{k) = g{r'^°^{Tk)) and f(fc) = ^(rfc), k = 
1, 2, 3. The subscript 2 ^ 3 at the end of the in- 
tegrand indicates that the contraction ^(2)^(3) 
6{t2 — Ts) is canceled by the identical contribu- 
tion coming from the average (C(''2)C('''3)) enter- 
ing the ((z/'"')^) in the last term of Eq. Re- 
peated application of the functional integration 
by part formula leads quickly to the following 
expression for the RHS of Eq. (jBip : 



/:^dr(l-e^)2 



+ ^(1 - e^fd^ ( 5(r<«)(r))(5(r<")(0) - r) 



The calculation follows the same lines of those 
leading to Eq. The averages can be calcu- 
lated using Eq. ^ for s(t) = r<°'(T) - r. It 
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is convenient to express the result in terms of 
Fourier components; for r = 0: 



-ibfc(l-e-r}exp{-- 



(-1 



As in the previous ID calculations we find a 

logarithmic divergence, this time at k 

and the leading contribution to the integral is 



2n 



dk 



{qkrf-liqkrr 



exp 



Using 50 / dxg{x) = 2r^ and ((t^'"')^) = q^/2, 
we obtain the final result, that is Eq. p3|) : 



(K')2)-i(i.2|o)(;' =-4(7re)-i/2inc5^ (B2) 
with c an 0(1) constant. 

B.2 The 3D case 

Calculation of (i^^ |0)c; in 3D is only slightly more 
involved. Using Eqs. pfTD)) and 1^, the RHS 
of Eq. (pij) can be expressed as the sum of three 
integrals: 

h = 3l r dri f" dT3 f " dT4 P 
x(^dpgc,^{l)do.g0^il)U^)U^)S{r{O) - r) 



3#4 

(B3) 



r dn f"^ dr2 r dT3 f " dr4 P 



x(l - e"i)e-"i+"^(%g„7(l)9a5/30(2 
x^^(l)C42K.(3)C.(4)5(r(0)-r)' 



3#4 



(B4) 



X 9^3507 (1)^7 (1)^43)?. (4)^(r(0) - r)] 



3#4 

(B5) 

where P = P({rfe}) = (1 - e"'i)e^-'+r4_ ^j^g 
ID case, the contractions between ^(3) and ^(4) 
are canceled by the last normalization term in 
Eq. jST]). 

In analogy with the analysis of Eq. ITSl) , the 
magnitude of the three integrals /2,2,3 can be 
estimated exploiting the fact that the integrals 
are dominated by ri,2 ^ e^^/^ <C 1 because 
of the factors (ti_2), and using everywhere 



.-1/2 



(then also P ~ r) and dr 



^^(3,4). The contractions with ^(1,2) lead to 
/ drfee'^''^(3, 4) 1. The same result is pro- 
duced from action of ,f (3, 4) on g{l)g{l,2)S{r'°\0) 
— r), that is evaluated with the functional in- 
tegration by part formula. We obtain in fact: 



1 In the same way we find: ^^(l) — *■ q{l—e^^)dj ~ 
qT2rv ~ 1, and ^^(2) ~ 1. Substituting into 
Eqs. (IB3IIB5P . we find /i,3 = 0{q^e-^^^) and 
/2 = 0(g2e-i). We thus obtain to Oiq'^e^^): 



z."|r = 0) = -^£r2dr(9„5/37(r'"'W) 

xa^g„^(r<«)(r))5(r<")(0))) 
= -T/^^'dr/d3r„ p[r^"^{T) = rjr<°)(0) = 0] 

X9a5/37(l'a)9/35a7(ra) 



Passing to Fourier components, and using Eq. 
l2l) we obtain the final result: 



xkH^{l - z^)z exp(-|k + 1| VV4), 

(B6) 

where = (k • l)/(fcZ), ((i/*'")^) = 3(jV2 and 
exp(-|k+l|2rV4) = Z,i/2,(rv(k+l), withZt(k) 
the generating function for r'°'(r) conditioned 
to r("'(0) = 0. 

The angular integral in Eq. (jB6p . is in the 

form /^^dz(l-z2)zexp(-(fc2+/2+2fcZz)T74) < 
0, and this tells us that the constant H in Eq. 
(p4|) is positive. 



Let us pass to the noise terms and consider first 
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